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Abstract. A major challenge in neuroscience is posed by the need for 
relating the emerging dynamical features of brain activity with the under¬ 
lying modular structure of neural connections, hierarchically organized 
throughout several scales. The spontaneous emergence of coherence and 
synchronization across such scales is crucial to neural function, while 
its anomalies often relate to pathological conditions. Here we provide 
a numerical study of synchronization dynamics in the human connec¬ 
tome network. Our purpose is to provide a detailed characterization of 
the recently uncovered broad dynamic regime, interposed between order 
and disorder, which stems from the hierarchical modular organization of 
the human connectome. In this regime -similar in essence to a Griffiths 
phase- synchronization dynamics are trapped within metastable attrac¬ 
tors of local coherence. Here we explore the role of noise, as an effective 
description of external perturbations, and discuss how its presence ac¬ 
counts for the ability of the system to escape intermittently from such 
attractors and explore complex dynamic repertoires of locally coherent 
states, in analogy with experimentally recorded patterns of cerebral ac¬ 
tivity. 
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1 Introduction 

The current mapping of neural connectivity patterns relies on advanced neuro¬ 
imaging techniques, which have recently allowed for the reconstruction of struc¬ 
tural human brain networks, establishing at an individual-based level which brain 
regions are mutually connected, as well as the strength of pairwise connections. 
The resulting “human connectome network” HE! has been found to be struc¬ 
tured in moduli or compartments -characterized by a much larger intra than 
inter connectivity- organized in a hierarchical fractal-like fashion across diverse 
scales mwm- On the other hand, functional connections between differ¬ 
ent brain regions can be inferred e.g. from correlations in neural activity as de¬ 
tected in EEG or fMRI time series. Unveiling the relation between structural and 




functional networks is a current challenge in modern neuroscience. In this con¬ 
text, a few pioneering works found that the hierarchical-modular organization 
of structural brain networks has remarkable implications for neural dynamics 
miojii]- As opposed to the case of simpler network structures, neural activ¬ 
ity propagates in hierarchical networks in a peculiar way. For example, models 
of neural activity propagation usually exhibit two familiar phases -percolating 
and non-percolating, respectively- but it has been recently found [12] that when 
such models operate on top of the “human connectome” structural network a 
novel intermediate regime, named a “Griffiths phase” mm emerges. This novel 
phase originates from the highly-diverse and relatively isolated structural mod¬ 
uli where dynamical activity may remain mostly localized for long time periods 

PP . 

Given that the correct brain functioning requires coherent neural activity 
at a wide range of scales mm, the study of synchronization among neural 
populations is one of the central ideas in computational neuroscience [ 17118 ] . 

In a recent work m, some of us scrutinized the special features of synchro¬ 
nization dynamics [20j using the canonical Kuramoto model for phase synchro¬ 
nization [211221281 . in the actual human connectome (HC) network [112124) . In 
analogy to what described above for activity propagation, we uncovered the ex¬ 
istence of a novel intermediate phase for synchronization dynamics, stemming 
from the hierarchical modular organization of the HC. Furthermore, we found 
that the dynamics in such a region presented a plethora of complex and inter¬ 
esting dynamical features [19] . 

Our goal here is to describe in more detail the complex behavior within such 
an intermediate regime, both in individual moduli and at a global brain level. 
We measure the fluctuations of the global order parameter as a function of the 
overall coupling strength, and we show that there is a broad region (rather than a 
unique “critical” point) with huge variability and response. Finally, we assess the 
role of noise and perturbations in the robustness of the metastable stated arising 
in the intermediate regime, and we show that adding intrinsic fluctuations to the 
picture of synchronization dynamics in hierarchical modular networks accounts 
for the ability of the brain to explore different attractors, giving access to the 
varied functional configurations recorded in experiments [ 25126127 ) . 

2 Kuramoto model in the Human-Connectome network 

The HC network we employ consists of a set of N = 998 nodes, each of them 
representing a population of neurons producing self-sustained oscillations [28], 
connected pairwise through a precise pattern of symmetric weighted edges, al¬ 
together determining a connectivity matrix W [Hi- 

On top of such a HC network, we implement a noisy Kuramoto dynamics, 
defined by the set of differential equations 121 [2212:11 : 

N 
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where 9i(t) is the phase at node i at time t, the intrinsic frequencies uj t - 
accounting for region heterogeneity- are extracted from some probability dis¬ 
tribution function g(w), Wij are the elements of the N x N weighted connec¬ 
tivity matrix W, k is an overall coupling parameter and r]i(t) is a zero-mean 
delta-correlated Gaussian noise, tuned by the real-valued amplitude a. 

The Kuramoto complex order parameter is defined by Z(t) = R(t)e = 
( e »0fc(t)) fe) where 0 < R < 1 gauges the overall coherence and ijj(t) is the average 
phase. It is common wisdom that for an (infinitely) large population of oscillators 
interacting in a fully connected network, the model exhibits a phase transition at 
some value of k, separating a coherent steady state (R > 0) from an incoherent 
one (R = 0, plus 1/y/N finite-size corrections) [21122123] . On the other hand, 
in the absence of frequency heterogeneity the system always reaches a coherent 
state. Thus, frequency heterogeneity is responsible for frustrating synchroniza¬ 
tion if the coupling strength is weak. Similarly, in our recent work |19j we argued 
that the combined effect of frequency heterogeneity and network heterogeneity 
(in particular, a hierarchical modular structure) can lead to much richer and 
interesting ways of ordering frustration. Here we explore that phenomenology in 
much deeper detail, introducing external stochastic fluctuations (i.e. noise) as 
the mechanism accounting for the ability of the system to explore metastable 
configurations. 


3 Results 

We considered the HC network |ll2] and employed standard community detec¬ 
tion algorithms [29|8| to identify the underlying modular structure. The optimal 
partition into communities i.e. the one maximizing the modularity parameter 
[50]— turns out to correspond to a division in 12 moduli [T9]. At a higher hier¬ 
archical level, a separation into just 2 moduli (roughly corresponding to the 2 
cerebral hemispheres) also provides a quite high modularity value. As illustrated 
in Fig. 1, 4 (out of the 12) moduli belong to one of the two hemispheres, 5 to 
the other, while 3 moduli (cyan, blue and red) overlap with both hemispheres. 
We label these two hierarchical levels as l = 2 (2 large moduli) and l = 1 (12 
smaller moduli), respectively. 

We have conducted computational analyses of the noisy Kuramoto model on 
top of the HC network and performed a number of new computational experi¬ 
ments complementing the analyses in our previous work |}19] . 

As illustrated in Fig. 2, beside the aforementioned coherent and incoherent 
phases (usually encountered in synchronizing systems) there is an intermediate 
regime between them exhibiting a large variability. Individual trajectories are 
depicted in the inset, for different values of the coupling strength k; observe in 
particular the irregular oscillations obtained for intermediate values of k. 

The reported values of ( R) t in the main plot of Fig. 2 correspond to the time- 
averaged value for a single realization in its steady state, considering up to a fixed 
maximum time T. The observed variability in the central region means either 
that (i) larger time windows would be required for the system to self-average 
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Fig. 1 . Adjacency matrix of the HC network with nodes ordered to emphasize its 
modular structure as highlighted by a community detection algorithm (see main text), 
showing also the partition into the 2 hemispheres (dashed lines). 12 moduli can be 
distinguished (each plotted with a different color); 4 of them correspond to one of the 
two hemispheres, 5 to the other, and only 3 moduli overlap with both hemispheres 
(cyan, blue and red moduli). Inter-modular connections (grey) are limited to small 
subsets, acting as interfaces or connectors between moduli. 


or (ii) that ergodicity is broken and for each parameter value the realization 
ends up in a different type of (stable or metastable) steady state, depending on 
the initial condition. This last possibility implies that the system may remain 
trapped in some sort of metastable states, from which it can escape away only 
after very rare and large fluctuations. 

These observations are robust against changes in the frequency distribution, 
connectivity matrix normalization, and other details, whereas the location and 
width of the intermediate phase are not universal. For example, Fig. 2 has been 
obtained for a Gaussian frequency distribution but similar curves are obtained 
for, usually employed, Lorentzian or uniform distributions. 

As this robust intermediate regime is reminiscent of Griffiths phases in net¬ 
works -posed in between order and disorder and emerging from rare-region ef¬ 
fects mm- it is natural to wonder how the structural network modularity 
affects synchronization dynamics in general. As a matter of fact, it is straightfor¬ 
ward to convince oneself that any network consisting of perfectly isolated moduli, 
each of them synchronized at different intrinsic frequencies and phases, should 
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control parameter k 


Fig. 2. (Main) Time-averaged value of the order parameter for the noisy Kuramoto 
dynamics running upon the Human Connectome network (998 nodes) with a Gaussian 
distribution of frequencies. Three different regimes emerge: an incoherent phase (k < 
1.6), a synchronous one (k > 3.3), and an intermediate irregular one. In this last, much 
larger averaging times would be required to obtain reliable mean values and these 
would depend upon initial conditions, reflecting metastability. We have chosen not so 
large measuring times (t = 100 for all values of fc) to illustrate the variability in the 
intermediate region. (Inset) Time-series for 4 different k values, indicated by arrows in 
the x-axis (from left to right: k = 0.5, 2.7, 3.0 and 4.0). 


exhibit oscillations of the collective order parameter, R, and these oscillations 
are preserved when the moduli are weakly interconnected m- Thus, in large 
networks without delays or other additional ingredients, time oscillations in the 
global coherence are the trademark of an underlying modular structure. 

To illustrate the role played by internal network modularity on global syn¬ 
chronization, Fig. 3 portraits the trajectories of the parameter Z(t) in the com¬ 
plex plane for different values of the control parameter fc, measured at different 
hierarchical levels: two (out of the existing 12) different small moduli (violet 
and orange curves), the two hemispheres (red and green), and the overall brain 
(blue). In the incoherent phase (panel a), the real and imaginary parts of Z 
fluctuate around zero at all scales in the hierarchy. On the other hand, in the 
coherent phase (panel d), all nodes are synchronized, and trajectories are circles 
with radii close to unity at all hierarchical levels 

A much richer behavior is found in the intermediate region: panel b (left) il¬ 
lustrates a situation in which one modulus (orange) is mostly coherent, while the 
other (violet) is not; however, hemispheres and global dynamics remain mostly 
unsynchronized (panel b (right)). In panel c (left), we have slightly increased 
the control parameter with respect to panel b, with a subsequent increase of 
the coherence for all hierarchical levels. Interestingly, as not all moduli exhibit 
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the same state of coherence, chaotic-like oscillations of the order parameter are 
observed at the global scale. 



Fig. 3. Phase portraits of the complex order parameter Z(t), measured at different 
scales in the hierarchy for a Gaussian g(w) (different realizations from those in Fig. 
2): two of the existing moduli are plotted in violet and orange, respectively, the two 
hemispheres in red and green, and the global scale in blue. Panels (a)-(d) correspond 
to values of the control parameter k = 1,3, 5 and 8, respectively (panels (b) and (c) 
have been split into two to enhance clarity), (a) In the non-ordered phase, the real and 
imaginary components of Z fluctuate around zero, not exhibiting synchronization at 
any scale, (b) In the early region of the intermediate phase, a few moduli are coherent 
(as the one in orange) but most of them remain unsynchronized (violet), and the 
system does not present coherence for upper scales in the hierarchy, (c) Increasing k, 
more heterogeneity of synchronization among moduli is found, and the system exhibits 
complex trajectories for the intermediate (hemispheres) and global scale, (d) In the 
coherent phase, all moduli are synchronized, and trajectories are concentric circles. 


We are interested in quantifying the observed variability of R in the interme¬ 
diate phase. To this end, we take a particular realization of frequencies (extracted 
from a Gaussian g(w)) and, starting from an initial -uniformly distributed- ran¬ 
dom configuration of individual phases, {9i(t = 0)}^, we measure the temporal 
standard deviation of the global coherence parameter R (after the transient) up 
to a maximum time T = 10 4 , 


a=(((R~(R) t ) 2 ) t y /2 


( 2 ) 
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as a function of the coupling strength k. Q 

As ergodicity may be broken, different initial conditions may lead to differ¬ 
ent attractors of the dynamics, therefore we also average a over 100 different 
independent realizations of the dynamical process. Results are illustrated in Fig. 
4, in which we also have plotted the diagram of the order parameter obtained 
for this particular realization of g(ui) averaged over the 100 realizations. Let us 
stress the following salient aspects: i) averaged time variabilities are small in the 
non-coherent (k < 1) as well as in the coherent (k > 5) phases, whereas much 
larger variabilities are found in the intermediate region (1 1$ k < 5); ii) the curve 
of time variabilities presents several peaks for the intermediate region, lying in 
the vicinity of values of the control parameter at which the system experiences 
a change in its level of coherence (see the corresponding jumps in the derivative 
of the order parameter); and finally, in) error bars are also larger in the inter¬ 
mediate phase; this variability of time variabilities means that different initial 
conditions can lead to different types of time-series, suggesting a large degree of 
metastability in the intermediate regime. 



control parameter k 


Fig. 4. Time-averaged order parameter R = (R(t)) and standard deviation of time- 
series, averaged over realizations with different -uniformly distributed- initial con¬ 
ditions. Maximal variability is found in the intermediate phase, where the system is 
neither too unsynchronized nor too coherent. Several peaks in the variability can be dis¬ 
tinguished (dashed lines), which appear at values of the control parameter k for which 
the system experiments a fast increase in global synchronization. Statistical sampling 
of different realizations indicate that errorbars are larger in the intermediate region, 
suggesting the existence of several attractors depending on the initial conditions. We 
have averaged 100 different realizations, each one integrated for 10 4 time steps. 


1 Notice that this definition of a, that we call, “time variability” is closely related to 
the chimera index introduced by Shanahan While chimera indices are averaged 
between individual network moduli and measure the onset of local coherence, a is 
defined at the global level and records fluctuations of the global order parameter. 
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3.1 Metastability in HMNs 


Our previous results vividly illustrate the existence of an intermediate region in 
which the HC exhibits maximal dynamical variability at the global scale, suggest¬ 
ing metastable behavior. In order to explore more directly whether metastable 
states exist, we now assess if the dynamics may present different attractors and, 
for some values of the control parameter k and noise amplitudes, if the system 
may switch between different global attractors with different levels of coherence. 

Fig. 5 (a) shows a time series of the global parameter, for a fixed realiza¬ 
tion of internal frequencies and initial phases. It clearly illustrates how the HC 
spontaneously switches between two different attractors. These type of events, 
however, are not easy to observe in the HC network. Due to the coarse-grained 
nature of the HC mapping, different attractors may actually have comparable 
average values of the coherence R , which makes their discrimination especially 
difficult at the global scale. 

Instead, such events are easier to spot in synthetic hierarchical modular net¬ 
works (HMN), such as proposed to model brain networks in an efficient way 
(see m and references therein). In such hierarchical networks, the effects of 
modularity and hierarchy are much enhanced, as they develop across a larger 
number of hierarchical levels than the one allowed by current imaging techniques 
for empirically obtained connectomes. 

All the previously reported phenomenology is still present in such HMNs (see 
j!9j): in particular, the phase diagram of the synchronization order parameter 
exhibits a phase transition with an intermediate region, where variability is much 
enhanced nn. Fig. 5 (b) illustrates the bi-stable nature of the global parameter 
in the intermediate phase for a HMN, in which metastability can be very well 
appreciated. This switching behavior closely resembles “up and down” states, 
which are well known to appear in certain phases of sleep or under anaesthesia 
(see [32] and refs, therein). 

We hypothesize that hierarchical modular networks in general (and the HC in 
particular) enable the possibility of a large repertoire of attractors, with different 
degrees of coherence and stability. Such metastability can be made evident and 
quantified by performing the following type of numerical test. Starting from a 
fixed random initial condition and considering a vanishing noise amplitude (i.e. 
a = 0), the system might deterministically fall into a number of different attrac¬ 
tors, each of them with an associated value of the global coherence depending on 
the initial conditions, the network structure, and the choice of natural frequen¬ 
cies. Once this attractor A is reached, the system is perturbed by switching on a 
non-vanishing noise amplitude (a > 0) during a finite time window. The system 
may remain stable in the same attractor A if the noise is weak enough (a -C 1). 
However, if larger values of the noise amplitude are chosen, the system may jump 
into another close, more stable, attractor. If the noise amplitude is very large 
(a 1), the system can in principle jump to any attractor, but, very likely, will 
also escape from it, wandering around a large fraction of the configuration space. 
After the perturbation time-window is over, we let the system relax once again, 
and check if the new resulting steady steady state B has changed with respect 


HMN 


(b) 

0.4 
0.3 
0.2 

0.4 
0.3 
0.2 

0.4 
0.3 
0.2 

0 20000 40000 60000 80000 

time 

Fig. 5. Time series exhibit metastability of the global synchronization in the Human 
Connectome and in HMNs, in the intermediate region, (a) Times series of the noisy 
Kuramoto dynamics in the HC with Gaussian g(w) in the intermediate region: for low 
noise amplitudes (a = 0 and 0.1), the system stays stable in the same attractor. But, for 
sufficiently large values, such as a = 0.2, the system is able to “jump” to another more 
coherent attractor, where it settles, (b) In HMNs (size N = 1024, 9 hierarchical levels), 
we observe the same phenomenology, but much enhanced: when noise is very low (a < 
0.45), the system tends to remain stable in a certain attractor (with a few exceptions 
after very large waiting times). Choosing a higher a (a > 0.5), the system exhibits 
bi-stable behavior, switching intermittently between two different attractors. For large 
enough a (a > 0.55), the dynamics becomes too erratic to appreciate metastability. 
Here, frequencies were extracted from a Lorentz distribution. 
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to A. In that case, we can conclude that the systems was in a metastable state 
A before the perturbation, and has reached another state B after it - potentially 
a metastable state itself. 

We have carried out this type of test using an artificial HMN (see Fig. 6) 
for a specific value of the control parameter, fc, belonging in the intermediate 
region. Natural frequencies are sampled from the a Lorentzian distribution g(uj) 
(as above, our main results are not sensible to this choice). Starting from a 
random initial configuration of phases, we integrate Eq. [l] up to time 500 with 
a = 0. After this, we introduce the external perturbation by switching the noise 
coefficient a to a certain non-zero value during a time window of duration 100. 
Finally we revert to a = 0 and continue the integration up to time t = 1000. 
The last steady state value is averaged over 10 4 realizations of initial conditions, 
networks, and intrinsic frequencies. 

As illustrated in Fig. 6, for low as well as for high values of the noise ampli¬ 
tude, the system has the same average order parameter close to {R)t,r uns — 0.2, 
as could have been anticipated. However, a resonant peak emerges for inter¬ 
mediate values of the noise, where the system switches to states with different 
levels of coherence. This plot explicitly illustrates the existence of metastability 
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and noise-induced jumps between attractors. As noise is enhanced, progressively 
more stable states are found, but above some noise threshold, the system does 
not remain trapped in a single attractor but jumps among many, resulting in a 
progressive decrease of the overall coherence. 


(a) 


(b) 



time 



Fig. 6. Perturbations can lead the system to more coherent attractors in the intermedi¬ 
ate non-coherent phase, (a) Order parameter R averaged in time over 10 4 realizations. 
A noise pulse of amplitude a is applied during the green interval. This same protocol 
is repeated for different values of a. (b) Average order parameter in the final steady 
state (after the noise pulse) as a function of a. For intermediate values of a, a resonant 
peak emerges for 1 < a < 10, illustrating that the system can jump to a close, more 
coherent on-average attractor. Simulations are run on HMN networks of size N = 1024, 
with 9 hierarchical levels. 


4 Discussion 

It is well established that in the absence of frequency dispersion, the Kuramoto 
dynamics leads to a perfectly coherent state, which is progressively achieved in 
time by following a bottom-up ordering dynamics in which increasingly larger 
communities become synchronized |33! . 

If a hierarchical modular networks is loosely connected, this type of “matryovska- 
doll” synchronization process is constrained at all levels by structural bottle¬ 
necks, bringing about anomalously-slow synchronization dynamics as recently 
reported PS- 

In the presence of intrinsic frequency dispersion the above slow ordering 
process is further frustrated m ■ For small values of k the system may re¬ 
main trapped into metastable states in which the loose connectivity between 
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some moduli does not allow them to overcome intrinsic-frequency differences 
and achieve coherence. While persistence in metastable states may extend indef¬ 
initely, experimental evidence suggests that the brain is able to switch between a 
rich repertoire of attractors We have shown that a simple description 

of neural coherence dynamics based on the noisy Kuramoto model may suf¬ 
fice to reproduce a very rich phenomenology, in hierarchical modular networks 
and in particular in the human connectome. The introduction of small fluctu¬ 
ations (exemplifying external perturbations, stimuli, or intrinsic stochasticity) 
allow the system to escape from metastable states and sample the configuration 
space, proving a paradigmatic modeling tool for the attractor surfing behavior 
suggested by experiments. Additional ingredients, such as explicit phase frus¬ 
tration [31j or time delays [28134] . should only add complexity to the structural 
frustration effect reported here, providing a finer description of brain activity. 
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